(defun is-square(x)
;  (declare (optimize (speed 3) (safety 0)))
 ;(declare (fixnum x))
  (let ((root (isqrt x)))
    (if (= x (expt root 2)) root nil)))

(defun find-solution(d)
;  (declare (optimize (speed 3) (safety 0)))
;  (declare (fixnum d))
  ;(if (is-square d) 0
  (loop for y from 1 to 10000
     for x = (is-square (1+ (* d (expt y 2))))
     when x do (return x)
     finally (return 0)))

(defun fractions(s)
  (loop 
;     for i from 1 to 10
     for m = 0 then (- (* d a) m)
     for d = 1 then (/ (- s (* m m)) d)
     for a = (isqrt s) then (floor (/ (+ (sqrt s) m) d))
     until (find (list m d a) values :test #'equal)
     for values = nil then (cons (list m d a) values)
     when (not (null values)) collect a ))

(defun finda(x)
;  (let* ((fracts (fractions x))
	; (terms (if (zerop (mod (length fracts) 2)) (reverse fracts) (reverse fracts))))
  (loop for h in (cddr (reverse (cons (isqrt x) (fractions x))))
     for n = 1 then d
     for d = n then x
     for x = (+ n (* h d)) then (+ (* h d) n)
     do (print (list h n d x))
     finally (return x)))


#|
( 97 106 109 139 149 151 157 163 166 172 181 193 199 202 211 214 233 241 244
 250 253 261 262 265 268 271 277 281 283 284 298 301 307 309 313 317 331 334
 337 343 349 353 358 365 367 373 379 382 388 393 394 397 406 409 412 417 419
 421 424 431 433 436 445 446 449 451 454 457 461 463 466 477 478 481 487 489
 491 493 501 502 508 509 511 517 521 523 524 526 533 537 538 541 547 549 550
 553 554 556 559 562 565 566 569 571 581 586 589 593 596 597 599 601 604 606
 607 610 613 614 617 619 622 628 629 631 633 634 637 641 643 647 649 652 653
 655 661 662 664 666 667 669 673 679 681 683 685 686 691 694 698 701 709 716
 718 719 721 724 733 737 739 742 746 749 751 753 754 757 758 763 764 766 769
 771 772 773 778 781 787 789 790 794 796 797 801 802 805 807 809 811 814 821
 823 826 829 834 835 838 844 845 849 853 856 857 859 861 862 863 865 869 871
 873 877 879 881 883 886 889 893 907 911 913 914 917 919 921 922 926 928 929
 931 932 937 941 946 947 949 951 953 954 955 956 958 964 965 967 970 971 972
 974 976 977 981 988 989 991 997 998 999 1000)
 
1766319049


public class test {

	
	static long doProcess(long x)
	{
		//;var x=document.ContFrac.firstNum.value;
		long xx=x;
		long xa=(long) java.lang.Math.sqrt(x);
		long a= xa;//parselong(xa+" ")
		if((a*a)==x)
		{
			return 0;
		}
		else
		{
			long[] num=new long[1000];
			num[0]=a;
			long b=a;
			long d=1;
			String s=a+";";
			String ss="";
			long j=0;
			long c=0;
			long r=0;
			long z=0;
			long ii=0;
			do
			{
				d=(long) Math.floor((x-(b*b))/d);
				c=(long) Math.floor((a+b)/d);
				r=(a+b)%d;
				b=Math.abs(r-a);
				s=s+c+",";
				num[(int) (j+1)]=c;
				j++;
			} while(c!=(2*a));
			long l=s.length();
			s=s.substring(0,(int) (l-1));
			//document.ContFrac.resultNum.value=s+"...("+j+" Repeating Digits)"
			if((j%2)==0)	//is even
			{
				ii=j-1;
			}
			else
			{
				ii=(2*j)-1;
				for(long m=1;m<=j;m++)
				{
					num[(int) (j+m)]=num[(int) m];
				}
			}
			//partial cinvergents
			ss="";
			long n=0;
			long h=0;
			x=0;
			
			n=1;
			d=num[(int) ii];
			for(long i=(ii-1);i>-1;i--)
			{
				h=num[(int) i];
				x=h*d+n;
				n=d;
				d=x;
				System.out.println("->h=" + h  + " n=" + n + " d=" + d + " x=" + x);
				
			}
			/*document.ContFrac.partial.value="("+d+")("+d+")"+"-"+"("+xx+")("+n+")("+n+") = 1"
			ss=ss+d+"/"+n+","
			l=ss.length
			ss=ss.substring(0,(l-1))
			document.ContFrac.equation.value=ss*/
			return d;
		}
	}	
	
	public static boolean isSquare(long x){
		double sqrt = Math.sqrt(x);
		return (sqrt == (double) (long) sqrt);
	}
	
	/**
	 * @param args
	 */
	public static void main(String[] args) {
//		long max = 1;
//		int imax = 1;
//		for (long i = 2; i<= 1000; i++){
//			if (! isSquare(i)){
//				long res = doProcess(i);
//				if (max < res ){
//					max = res;
//					imax = (int) i;
//				}
//			}
//		}
//		System.out.println(imax);
//		System.out.println(java.lang.Long.MAX_VALUE);
		System.out.println(doProcess(86));
	}
}

|#

;;; Looked up Pell's equation on Wikipedia...
(defun approximate-root (x n &optional (a 1) (b 0))
  (if (= n 0)
      (truncate (+ (* a (sqrt x)) b))
      (let ((c (truncate (+ (* a (sqrt x)) b))))
	(let ((a (/ a (- (* a a x) (expt (- b c) 2))))
	      (b (/ (- c b) (- (* a a x) (expt (- b c) 2)))))
	  (+ c (/ (approximate-root x (1- n) a b)))))))

(defun solve-for (n)
  (do* ((i 0 (1+ i))
	(root (approximate-root n i) (approximate-root n i)))
       ((= (- (expt (numerator root) 2)
	      (* n (expt (denominator root) 2)))
	   1)
	(numerator root))))

(do ((i 1 (1+ i))
     (max 0) best)
    ((> i 1000) best)
  (when (not (is-square i))
    (let ((solution (solve-for i)))
      (when (> solution max)
	(setf max solution best i)))))
